Time Series Forecasting
Bianca Gunawan
2024-05-06
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(base)
library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, lastlibrary(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooData Preprocessing
dataPath <- "/Users/bianca/Documents/U of Chicago/Classes/Time Series"
contracts_volume <- read.csv(paste(dataPath,'Contracts_Volume.csv',sep = '/'), header=TRUE)
contracts_class <- read.csv(paste(dataPath,'Contracts_Classification.csv',sep = '/'), header=TRUE)Computing Floor Volume
contracts_volume$Total.Volume <- as.numeric(contracts_volume$Total.Volume)
contracts_volume$Electronic.Volume <- as.numeric(gsub(",", "", contracts_volume$Electronic.Volume))
contracts_volume$Floor.Volume <- contracts_volume$Total.Volume - contracts_volume$Electronic.VolumeJoin Contracts Class and Volume
joined_data <- contracts_volume %>%
left_join(contracts_class, by = c("Commodity.Indicator" = "Commodity.Code"))Sort Volume Data Based on Division
cme_data <- joined_data
cme_data$Division <- "CME"
cme_data$Date <- as.Date(cme_data$Date, format= "%m/%d/%Y")
imm_data <- filter(joined_data, Division == "IMM")
imm_data$Date <- as.Date(imm_data$Date, format= "%m/%d/%Y")
iom_data <- filter(joined_data, Division == "IOM")
iom_data$Date <- as.Date(iom_data$Date, format= "%m/%d/%Y")Aggregate Volume by Commodity Code and Month
cme_agg <- cme_data %>%
group_by(Commodity.Indicator, Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
imm_agg <- imm_data %>%
group_by(Commodity.Indicator, Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
iom_agg <- iom_data %>%
group_by(Commodity.Indicator, Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
print(cme_agg)## # A tibble: 21,048 × 5
## # Groups: Commodity.Indicator [559]
## Commodity.Indicator Date Total_Electronic.Volume Total_Total.Volume
## <chr> <date> <dbl> <dbl>
## 1 16E 2011-02-01 11 11
## 2 16E 2011-03-01 24 24
## 3 16E 2011-04-01 85 85
## 4 16E 2011-05-01 63 63
## 5 16E 2011-06-01 70 70
## 6 16E 2011-07-01 44 44
## 7 16E 2011-08-01 8 8
## 8 16E 2011-09-01 12 12
## 9 16E 2011-10-01 5 5
## 10 1A 2007-05-01 0 116
## # ℹ 21,038 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>print(imm_agg)## # A tibble: 5,723 × 5
## # Groups: Commodity.Indicator [78]
## Commodity.Indicator Date Total_Electronic.Volume Total_Total.Volume
## <chr> <date> <dbl> <dbl>
## 1 16E 2011-02-01 11 11
## 2 16E 2011-03-01 24 24
## 3 16E 2011-04-01 85 85
## 4 16E 2011-05-01 63 63
## 5 16E 2011-06-01 70 70
## 6 16E 2011-07-01 44 44
## 7 16E 2011-08-01 8 8
## 8 16E 2011-09-01 12 12
## 9 16E 2011-10-01 5 5
## 10 36E 2011-03-01 27 27
## # ℹ 5,713 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>print(iom_agg)## # A tibble: 16,399 × 5
## # Groups: Commodity.Indicator [475]
## Commodity.Indicator Date Total_Electronic.Volume Total_Total.Volume
## <chr> <date> <dbl> <dbl>
## 1 1A 2007-05-01 0 116
## 2 1A 2008-01-01 0 25
## 3 1A 2009-04-01 0 14
## 4 1A 2011-03-01 76 426
## 5 1A 2011-06-01 26 26
## 6 1A 2011-07-01 0 220
## 7 1A 2011-08-01 11 11
## 8 1A 2011-09-01 4 4
## 9 1A 2011-11-01 146 146
## 10 1A 2011-12-01 25 28
## # ℹ 16,389 more rows
## # ℹ 1 more variable: Total_Floor.Volume <dbl>Aggregate by Month only
cme_agg2 <- cme_data %>%
group_by(Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
imm_agg2 <- imm_data %>%
group_by(Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))
iom_agg2 <- iom_data %>%
group_by(Date) %>%
summarise(Total_Electronic.Volume = sum(Electronic.Volume),
Total_Total.Volume = sum(Total.Volume),
Total_Floor.Volume = sum(Floor.Volume))Loading Price Data
imputed_cme <- read.csv(paste(dataPath,'imputed_cme.csv',sep = '/'), header=TRUE)
names(imputed_cme) <- c("Date", "Price")
imputed_cme$Date <- as.Date(imputed_cme$Date)
imputed_imm <- read.csv(paste(dataPath,'imputed_imm.csv',sep = '/'), header=TRUE)
names(imputed_imm) <- c("Date", "Price")
imputed_imm$Date <- as.Date(imputed_imm$Date)
imputed_iom <- read.csv(paste(dataPath,'imputed_iom.csv',sep = '/'), header=TRUE)
names(imputed_iom) <- c("Date", "Price")
imputed_iom$Date <- as.Date(imputed_iom$Date)EDA
cme_summary <- left_join(imputed_cme, cme_agg2, by = "Date")
imm_summary <- left_join(imputed_imm, imm_agg2, by = "Date")
iom_summary <- left_join(imputed_iom, iom_agg2, by = "Date")
cme_numeric_data <- cme_summary %>%
select_if(is.numeric)
cme_cor_matrix <- cor(cme_numeric_data, use = "complete.obs", method = "pearson")
cme_cor <- as.data.frame(cme_cor_matrix)
imm_numeric_data <- imm_summary %>%
select_if(is.numeric)
imm_cor_matrix <- cor(imm_numeric_data, use = "complete.obs", method = "pearson")
imm_cor <- as.data.frame(imm_cor_matrix)
iom_numeric_data <- iom_summary %>%
select_if(is.numeric)
iom_cor_matrix <- cor(iom_numeric_data, use = "complete.obs", method = "pearson")
iom_cor <- as.data.frame(iom_cor_matrix)
print(cme_cor)## Price Total_Electronic.Volume Total_Total.Volume
## Price 1.00000000 0.6530748 0.6802748
## Total_Electronic.Volume 0.65307484 1.0000000 0.9736634
## Total_Total.Volume 0.68027480 0.9736634 1.0000000
## Total_Floor.Volume -0.06265014 -0.3726483 -0.1512653
## Total_Floor.Volume
## Price -0.06265014
## Total_Electronic.Volume -0.37264828
## Total_Total.Volume -0.15126533
## Total_Floor.Volume 1.00000000print(imm_cor)## Price Total_Electronic.Volume Total_Total.Volume
## Price 1.0000000 0.3287382 0.4764459
## Total_Electronic.Volume 0.3287382 1.0000000 0.9465020
## Total_Total.Volume 0.4764459 0.9465020 1.0000000
## Total_Floor.Volume 0.2972396 -0.4706299 -0.1607259
## Total_Floor.Volume
## Price 0.2972396
## Total_Electronic.Volume -0.4706299
## Total_Total.Volume -0.1607259
## Total_Floor.Volume 1.0000000print(iom_cor)## Price Total_Electronic.Volume Total_Total.Volume
## Price 1.00000000 -0.02515719 0.05596529
## Total_Electronic.Volume -0.02515719 1.00000000 0.98393410
## Total_Total.Volume 0.05596529 0.98393410 1.00000000
## Total_Floor.Volume 0.43943993 -0.29485138 -0.11951925
## Total_Floor.Volume
## Price 0.4394399
## Total_Electronic.Volume -0.2948514
## Total_Total.Volume -0.1195193
## Total_Floor.Volume 1.0000000It looks like Floor Volume is not a good indicator for CME class seats prices. Both Electronic Volume and Total Volume are not good indicators for IOM class seats prices. We can also see that Electronic Volume and Total Volume are highly correlated in all seat classes, so only one of those will be used as predictor.
Task A: Algorithms
Split train and test
train_data_cme <- subset(cme_summary, Date < "2013-01-01")
test_data_cme <- subset(cme_summary, Date >= "2013-01-01")
train_data_imm <- subset(imm_summary, Date < "2013-01-01")
test_data_imm <- subset(imm_summary, Date >= "2013-01-01")
train_data_iom <- subset(iom_summary, Date < "2013-01-01")
test_data_iom <- subset(iom_summary, Date >= "2013-01-01")
price_cme <- ts(train_data_cme$Price, frequency = 12)
price_imm <- ts(train_data_imm$Price, frequency = 12)
price_iom <- ts(train_data_iom$Price, frequency = 12)Linear Regression
lm_cme <- lm(Price ~ Total_Total.Volume, data = train_data_cme)
lm_forecast_cme <- predict(lm_cme, newdata = test_data_cme, n.ahead = 12)
summary(lm_cme)##
## Call:
## lm(formula = Price ~ Total_Total.Volume, data = train_data_cme)
##
## Residuals:
## Min 1Q Median 3Q Max
## -419623 -112256 -7097 82349 791934
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.166e+05 4.127e+04 2.826 0.0054 **
## Total_Total.Volume 2.945e-03 2.371e-04 12.420 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 197600 on 142 degrees of freedom
## Multiple R-squared: 0.5207, Adjusted R-squared: 0.5173
## F-statistic: 154.3 on 1 and 142 DF, p-value: < 2.2e-16lm_imm <- lm(Price ~ Total_Total.Volume, data = train_data_imm)
lm_forecast_imm <- predict(lm_imm, newdata = test_data_imm, n.ahead = 12)
lm_iom <- lm(Price ~ Total_Floor.Volume, data = train_data_iom)
lm_forecast_iom <- predict(lm_iom, newdata = test_data_iom, n.ahead = 12)For CME linear regression, Electronic Volume and Total Volume are correlated so only Total Volume is used. Floor Volume doesn’t have high correltion with price so it’s not used as predictor. For IMM linear regression, Electronic Volume and Floor Volume is picked instead so there’s no overlap between Total Volume and Floor Volume as predictors. For IOM, only Floor Volume is highly correlated with price so it’s the only predictor.
Linear regression with ARMA errors (use arima with xreg)
arma_model_cme <- auto.arima(price_cme, xreg = train_data_cme$Total_Total.Volume)
arma_forecast_cme <- forecast(arma_model_cme, xreg = test_data_cme$Total_Total.Volume, h=12)
arma_model_imm <- auto.arima(price_imm, xreg = train_data_imm$Total_Total.Volume)
arma_forecast_imm <- forecast(arma_model_imm, xreg = test_data_imm$Total_Total.Volume, h=12)
arma_model_iom <- auto.arima(price_iom, xreg = train_data_iom$Total_Total.Volume)
arma_forecast_iom <- forecast(arma_model_iom, xreg = test_data_iom$Total_Total.Volume, h=12)Holt-Winters
library(forecast)
hw_model_cme <- hw(price_cme, xreg = train_data_cme$Total_Total.Volume)
hw_forecast_cme <- forecast(hw_model_cme, xreg = test_data_cme$Total_Total.Volume, h = 12)
hw_model_imm <- hw(price_imm,xreg = train_data_imm$Total_Total.Volume)
hw_forecast_imm <- forecast(hw_model_imm, xreg = test_data_imm$Total_Total.Volume, h = 12)
hw_model_iom <- hw(price_iom, xreg = train_data_iom$Total_Total.Volume)
hw_forecast_iom <- forecast(hw_model_iom, xreg = test_data_iom$Total_Total.Volume, h = 12)ARIMA
arima_cme <- auto.arima(price_cme, seasonal = FALSE)
arima_forecast_cme <- forecast(arima_cme, h = 12)
arima_imm <- auto.arima(price_imm, seasonal = FALSE)
arima_forecast_imm <- forecast(arima_imm, h = 12)
arima_iom <- auto.arima(price_iom, seasonal = FALSE)
arima_forecast_iom <- forecast(arima_iom, h = 12)Seasonal ARIMA
# price_cme, price_imm, price_iom already has frequency of 12
sarima_cme <- auto.arima(price_cme, seasonal = TRUE)
sarima_forecast_cme <- forecast(sarima_cme, h = 12)
sarima_imm <- auto.arima(price_imm, seasonal = TRUE)
sarima_forecast_imm <- forecast(sarima_imm, h = 12)
sarima_iom <- auto.arima(price_iom, seasonal = TRUE)
sarima_forecast_iom <- forecast(sarima_iom, h = 12)Fractional ARIMA
# Looking at ACF
acf(price_cme)
acf(price_imm)
acf(price_iom)
The autocorrelations show trend, so arfima is applicable.
library(fracdiff)
d_cme <- fracdiff(price_cme)
st_cme <- diffseries(price_cme, d_cme$d)
arfima_cme <- auto.arima(st_cme, xreg = train_data_cme$Total_Total.Volume)
arfima_forecast_cme <- forecast(arfima_cme, xreg = test_data_cme$Total_Total.Volume, h = 12)
d_imm <- fracdiff(price_imm)
st_imm <- diffseries(price_imm, d_imm$d)
arfima_imm <- auto.arima(st_imm, xreg = train_data_imm$Total_Total.Volume)
arfima_forecast_imm <- forecast(arfima_imm, xreg = test_data_imm$Total_Total.Volume, h = 12)
d_iom <- fracdiff(price_iom)
st_iom <- diffseries(price_iom, d_iom$d)
arfima_iom <- auto.arima(st_iom, xreg = train_data_iom$Total_Total.Volume)
arfima_forecast_iom <- forecast(arfima_iom, xreg = test_data_iom$Total_Total.Volume, h = 12)ARMA and GARCH
#install.packages("fGarch")
library(fGarch)## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")num_train_data_cme <- scale(select_if(train_data_cme, is.numeric))
num_train_data_imm <- scale(select_if(train_data_imm, is.numeric))
num_train_data_iom <- scale(select_if(train_data_iom, is.numeric))
num_test_data_cme <- scale(select_if(test_data_cme, is.numeric))
num_test_data_imm <- scale(select_if(test_data_imm, is.numeric))
num_test_data_iom <- scale(select_if(test_data_iom, is.numeric))
garch_model_cme <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_cme, trace = FALSE)## Warning in sqrt(diag(fit$cvar)): NaNs producedgarch_forecast_cme <- predict(garch_model_cme, data = num_test_data_cme, n.ahead = 12)
garch_model_imm <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_imm, trace = FALSE)
garch_forecast_imm <- predict(garch_model_imm, data = num_test_data_imm, n.ahead = 12)
garch_model_iom <- garchFit(Price ~ arma(1,1) + garch(1,1), data = num_train_data_iom, trace = FALSE)
garch_forecast_iom <- predict(garch_model_iom, data = num_test_data_iom, n.ahead = 12)Compare results sMAPE across models
cme_rmse <- data.frame(Commodity = "CME",
LM = accuracy(lm_forecast_cme, test_data_cme$Price)[[2]],
ARMA = accuracy(arma_forecast_cme, test_data_cme$Price)[[2]],
HW = accuracy(hw_forecast_cme, test_data_cme$Price)[[2]],
ARIMA = accuracy(arima_forecast_cme, test_data_cme$Price)[[2]],
SARIMA = accuracy(sarima_forecast_cme, test_data_cme$Price)[[2]],
ARFIMA = accuracy(arfima_forecast_cme, test_data_cme$Price)[[2]],
GARCH = accuracy(garch_forecast_cme$meanForecast, num_test_data_cme[,1])[[2]])
imm_rmse <- data.frame(Commodity = "IMM",
LM = accuracy(lm_forecast_imm, test_data_imm$Price)[[2]],
ARMA = accuracy(arma_forecast_imm, test_data_imm$Price)[[2]],
HW = accuracy(hw_forecast_imm, test_data_imm$Price)[[2]],
ARIMA = accuracy(arima_forecast_imm, test_data_imm$Price)[[2]],
SARIMA = accuracy(sarima_forecast_imm, test_data_imm$Price)[[2]],
ARFIMA = accuracy(arfima_forecast_imm, test_data_imm$Price)[[2]],
GARCH = accuracy(garch_forecast_imm$meanForecast, num_test_data_imm[,1])[[2]])
iom_rmse <- data.frame(Commodity = "IOM",
LM = accuracy(lm_forecast_iom, test_data_iom$Price)[[2]],
ARMA = accuracy(arma_forecast_iom, test_data_iom$Price)[[2]],
HW = accuracy(hw_forecast_iom, test_data_iom$Price)[[2]],
ARIMA = accuracy(arima_forecast_iom, test_data_iom$Price)[[2]],
SARIMA = accuracy(sarima_forecast_iom, test_data_iom$Price)[[2]],
ARFIMA = accuracy(arfima_forecast_iom, test_data_iom$Price)[[2]],
GARCH = accuracy(garch_forecast_iom$meanForecast, num_test_data_iom[,1])[[2]])
combined_rmse <- rbind(cme_rmse, imm_rmse, iom_rmse)
combined_rmse$Best_Model <- apply(combined_rmse[, -1], 1, function(x) {
colnames(combined_rmse[, -1])[which.min(abs(x))]
})
print(combined_rmse)## Commodity LM ARMA HW ARIMA SARIMA ARFIMA
## 1 CME 271928.7 -241806.2206 -192.5793 31291.667 31291.667 512777.13
## 2 IMM 242626.0 25469.1639 33720.7608 34281.914 34281.914 268904.23
## 3 IOM 116316.6 -900.4912 14852.5595 -1071.768 -4836.629 89559.59
## GARCH Best_Model
## 1 1.064806 GARCH
## 2 1.242325 GARCH
## 3 1.606774 GARCHBest Model: - GARCH for all
Task B: BSTS
Load data
sts <- read.csv(paste(dataPath,'sts_data_2023.csv',sep = '/'), header=TRUE)
sts <- na.omit(sts)
sts_1990 <- sts %>%
filter(as.numeric(format(as.Date(Date, "%m/%d/%y"), "%y")) > 23) %>%
mutate(Date = format(as.Date(Date, "%m/%d/%y"), "19%y-%m-%d"))
sts_2000 <- sts %>%
filter(as.numeric(format(as.Date(Date, "%m/%d/%y"), "%y")) <= 23) %>%
mutate(Date = format(as.Date(Date, "%m/%d/%y"), "20%y-%m-%d"))
sts <- rbind(sts_1990, sts_2000)
sts <- xts(sts[,-1], order.by = as.Date(sts$Date))
#install.packages("bsts")
library(bsts)## Loading required package: BoomSpikeSlab## Loading required package: Boom##
## Attaching package: 'Boom'## The following object is masked from 'package:stats':
##
## rWishart##
## Attaching package: 'BoomSpikeSlab'## The following object is masked from 'package:stats':
##
## knots##
## Attaching package: 'bsts'## The following object is masked from 'package:BoomSpikeSlab':
##
## SuggestBurndata_fit1 <- window(sts, end = "1979-06-30")
data_fit2 <- window(sts,end = "2004-06-30")
data_fit3 <- window(sts, end = "2022-12-31")
new_data1 <- head(window(sts, start = "1979-09-30"), 5)
new_data1_CL <- new_data1[, -3]
new_data1_GC <- new_data1[, -1]
new_data2 <- head(window(sts, start = "2004-09-30"), 5)
new_data2_CL <- new_data1[, -3]
new_data2_GC <- new_data1[, -1]
new_data3 <- head(window(sts, start = "2023-03-31"), 1)
new_data_3CL <- new_data1[, -3]
new_data3_GC <- new_data1[, -1] Fit Models for CL
Dynamic regression of other 3 variables (lags)
# FIT 1
ss <- AddLocalLevel(list(), data_fit1$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit1$Mean_CL ~ data_fit1$Mean_GC + data_fit1$Mean_CPI + data_fit1$Mean_SPX)
model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit1, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:43 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:43 2024 =-=-=-=-=plot(model_dynamic)
pred_CL_dynamic_1 <- predict.bsts(model_dynamic, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_dynamic_1)
# FIT 2
ss <- AddLocalLevel(list(), data_fit2$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit2$Mean_CL ~ data_fit2$Mean_GC + data_fit2$Mean_CPI + data_fit2$Mean_SPX)
model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit2, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:44 2024 =-=-=-=-=plot(model_dynamic)
pred_CL_dynamic_2 <- predict.bsts(model_dynamic, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_dynamic_2)
# FIT 3
ss <- AddLocalLevel(list(), data_fit3$Mean_CL)
ss <- AddDynamicRegression(ss, data_fit3$Mean_CL ~ data_fit3$Mean_GC + data_fit3$Mean_CPI + data_fit3$Mean_SPX)
model_dynamic <- bsts(Mean_CL ~ ., state.specification = ss, data = data_fit3, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:44 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=plot(model_dynamic)
pred_CL_dynamic_3 <- predict.bsts(model_dynamic, newdata = new_data3, horizon = 1, quantiles = c(.025, .975))
plot(pred_CL_dynamic_3)
Different trend component for each model
# FIT 1
ss <- AddLocalLinearTrend(list(), data_fit1$Mean_CL)
model_trend <- bsts(data_fit1$Mean_CL, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=plot(model_trend)
pred_CL_trend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_1)
# FIT 2
ss <- AddLocalLinearTrend(list(), data_fit2$Mean_CL)
model_trend <- bsts(data_fit2$Mean_CL, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=plot(model_trend)
pred_CL_trend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_2)
# FIT 3
ss <- AddLocalLinearTrend(list(), data_fit3$Mean_CL)
model_trend <- bsts(data_fit3$Mean_CL, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:45 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:45 2024 =-=-=-=-=plot(model_trend)
pred_CL_trend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_trend_3)
State-specification with semi-local linear trend
# FIT 1
ss <- AddSemilocalLinearTrend(list(), data_fit1$Mean_CL)
model_trend <- bsts(data_fit1$Mean_CL, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=plot(model_trend)
pred_CL_semitrend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_1)
# FIT 2
ss <- AddSemilocalLinearTrend(list(), data_fit2$Mean_CL)
model_trend <- bsts(data_fit2$Mean_CL, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=plot(model_trend)
pred_CL_semitrend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_2)
# FIT 3
ss <- AddSemilocalLinearTrend(list(), data_fit3$Mean_CL)
model_trend <- bsts(data_fit3$Mean_CL, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:46 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:46 2024 =-=-=-=-=plot(model_trend)
pred_CL_semitrend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_CL_semitrend_3)
Fit Models for GC
Dynamic regression of other 3 variables (lags)
# FIT 1
ss <- AddLocalLevel(list(), data_fit1$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit1$Mean_GC ~ data_fit1$Mean_CL + data_fit1$Mean_CPI + data_fit1$Mean_SPX)
model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit1, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:47 2024 =-=-=-=-=plot(model_dynamic)
pred_GC_dynamic_1 <- predict.bsts(model_dynamic, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_dynamic_1)
# FIT 2
ss <- AddLocalLevel(list(), data_fit2$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit2$Mean_GC ~ data_fit2$Mean_CL + data_fit2$Mean_CPI + data_fit2$Mean_SPX)
model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit2, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:47 2024 =-=-=-=-=plot(model_dynamic)
pred_GC_dynamic_2 <- predict.bsts(model_dynamic, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_dynamic_2)
# FIT 3
ss <- AddLocalLevel(list(), data_fit3$Mean_GC)
ss <- AddDynamicRegression(ss, data_fit3$Mean_GC ~ data_fit3$Mean_CL + data_fit3$Mean_CPI + data_fit3$Mean_SPX)
model_dynamic <- bsts(Mean_GC ~ ., state.specification = ss, data = data_fit3, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:47 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=plot(model_dynamic)
pred_GC_dynamic_3 <- predict.bsts(model_dynamic, newdata = new_data3, horizon = 1, quantiles = c(.025, .975))
plot(pred_GC_dynamic_3)
Different trend component for each model
# FIT 1
ss <- AddLocalLinearTrend(list(), data_fit1$Mean_GC)
model_trend <- bsts(data_fit1$Mean_GC, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=plot(model_trend)
pred_GC_trend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_1)
# FIT 2
ss <- AddLocalLinearTrend(list(), data_fit2$Mean_GC)
model_trend <- bsts(data_fit2$Mean_GC, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:48 2024 =-=-=-=-=plot(model_trend)
pred_GC_trend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_2)
# FIT 3
ss <- AddLocalLinearTrend(list(), data_fit3$Mean_GC)
model_trend <- bsts(data_fit3$Mean_GC, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:48 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=plot(model_trend)
pred_GC_trend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_trend_3)
State-specification with semi-local linear trend
# FIT 1
ss <- AddSemilocalLinearTrend(list(), data_fit1$Mean_GC)
model_trend <- bsts(data_fit1$Mean_GC, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=plot(model_trend)
pred_GC_semitrend_1 <- predict.bsts(model_trend, newdata = new_data1, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_1)
# FIT 2
ss <- AddSemilocalLinearTrend(list(), data_fit2$Mean_GC)
model_trend <- bsts(data_fit2$Mean_GC, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=plot(model_trend)
pred_GC_semitrend_2 <- predict.bsts(model_trend, newdata = new_data2, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_2)
# FIT 3
ss <- AddSemilocalLinearTrend(list(), data_fit3$Mean_GC)
model_trend <- bsts(data_fit3$Mean_GC, state.specification = ss, niter = 200)## =-=-=-=-= Iteration 0 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 20 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 40 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 60 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 80 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 100 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 120 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 140 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 160 Mon Oct 28 13:41:49 2024 =-=-=-=-=
## =-=-=-=-= Iteration 180 Mon Oct 28 13:41:49 2024 =-=-=-=-=plot(model_trend)
pred_GC_semitrend_3 <- predict.bsts(model_trend, newdata = new_data3, horizon = 5, quantiles = c(.025, .975))
plot(pred_GC_semitrend_3)
Testing RMSE
CL_fit1_rmse <- data.frame(Commodity = "CL Fit 1",
DynamicRegression = accuracy(pred_CL_dynamic_1$mean, new_data1$Mean_CL)[[2]],
LinearTrend = accuracy(pred_CL_trend_1$mean, new_data1$Mean_CL)[[2]],
SemiLocalTrend = accuracy(pred_CL_semitrend_1$mean, new_data1$Mean_CL)[[2]])
CL_fit2_rmse <- data.frame(Commodity = "CL Fit 2",
DynamicRegression = accuracy(pred_CL_dynamic_2$mean, new_data2$Mean_CL)[[2]],
LinearTrend = accuracy(pred_CL_trend_2$mean, new_data2$Mean_CL)[[2]],
SemiLocalTrend = accuracy(pred_CL_semitrend_1$mean, new_data2$Mean_CL)[[2]])
CL_fit3_rmse <- data.frame(Commodity = "CL Fit 3",
DynamicRegression = accuracy(pred_CL_dynamic_3$mean, new_data3$Mean_CL)[[2]],
LinearTrend = accuracy(pred_CL_trend_3$mean, new_data3$Mean_CL)[[2]],
SemiLocalTrend = accuracy(pred_CL_semitrend_3$mean, new_data3$Mean_CL)[[2]])
GC_fit1_rmse <- data.frame(Commodity = "GC Fit 1",
DynamicRegression = accuracy(pred_GC_dynamic_1$mean, new_data1$Mean_GC)[[2]],
LinearTrend = accuracy(pred_GC_trend_1$mean, new_data1$Mean_GC)[[2]],
SemiLocalTrend = accuracy(pred_GC_semitrend_1$mean, new_data1$Mean_GC)[[2]])
GC_fit2_rmse <- data.frame(Commodity = "GC Fit 2",
DynamicRegression = accuracy(pred_GC_dynamic_2$mean, new_data2$Mean_GC)[[2]],
LinearTrend = accuracy(pred_GC_trend_2$mean, new_data2$Mean_GC)[[2]],
SemiLocalTrend = accuracy(pred_GC_semitrend_1$mean, new_data2$Mean_GC)[[2]])
GC_fit3_rmse <- data.frame(Commodity = "GC Fit 3",
DynamicRegression = accuracy(pred_GC_dynamic_3$mean, new_data3$Mean_GC)[[2]],
LinearTrend = accuracy(pred_GC_trend_3$mean, new_data3$Mean_GC)[[2]],
SemiLocalTrend = accuracy(pred_GC_semitrend_3$mean, new_data3$Mean_GC)[[2]])
combined_rmse <- rbind(CL_fit1_rmse, CL_fit2_rmse, CL_fit3_rmse, GC_fit1_rmse, GC_fit2_rmse, GC_fit3_rmse)
combined_rmse$Best_Model <- apply(combined_rmse[, -1], 1, function(x) {
colnames(combined_rmse[, -1])[which.min(abs(x))]
})
print(combined_rmse)## Commodity DynamicRegression LinearTrend SemiLocalTrend Best_Model
## 1 CL Fit 1 30.39734 6.255817 12.663403 LinearTrend
## 2 CL Fit 2 38.35451 3.689569 28.469102 LinearTrend
## 3 CL Fit 3 78.29768 2.932757 5.081877 LinearTrend
## 4 GC Fit 1 458.76358 249.223872 266.072375 LinearTrend
## 5 GC Fit 2 347.42962 11.774606 134.081005 LinearTrend
## 6 GC Fit 3 1531.99188 195.287713 204.353104 LinearTrendQUESTIONS Local linear trend works best with everything except fit 3 of GC, which is best modeled using Semi Local Linear Trend state-specification. Local Linear Trend is best overall.
Lags are not useful in predicting CL and GC. All coefficients are 0 meaning no meaningful linear relationship between variables.